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ABSTRACT 

We investigate the process of metal-free star formation in the first galaxies with 
a high-resolution cosmological simulation. We consider the cosmologically motivated 
scenario in which a strong molecule-destroying Lyman- Werner (LW) background in- 
hibits effective cooling in low-mass haloes, delaying star formation until the collapse 
or more massive haloes. Only when molecular hydrogen (H2) can self-shield from LW 
radiation, which requires a halo capable of cooling by atomic line emission, will star 
formation be possible. To follow the formation of multiple gravitationally bound ob- 
jects, at high gas densities we introduce sink particles which accrete gas directly from 
the computational grid. We find that in a 1 Mpc 3 (comoving) box, runaway collapse 
first occurs in a 3 x 10 7 Mq dark matter halo at z ~ 12 assuming a background inten- 
sity of J 21 = 100. Due to a runaway increase in the H2 abundance and cooling rate, a 
self-shielding, supersonically turbulent core develops abruptly with ~ 10 4 Mq in cold 
gas available for star formation. We analyze the formation of this self-shielding core, 
the character of turbulence, and the prospects for star formation. Due to a lack of frag- 
mentation on scales we resolve, we argue that LW-delayed metal-free star formation 
in atomic cooling haloes is very similar to star formation in primordial minihaloes, 
although in making this conclusion we ignore internal stellar feedback. Finally, we 
briefly discuss the detectability of metal-free stellar clusters with the James Webb 
Space Telescope. 
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1 INTRODUCTION 

Theoretical studies have suggested that the first stars, Pop- 
ulation III (Pop III), formed in ~ 10 6 Mq dark matter 'mini- 
haloes' at redshifts z ~ 15 — 40 (Couchman & Rees 1986 



Haiman et al.| [T996 Teg mark et al.||1997[ ). Due to a lack of 



efficient coolants in metal-free gas, Pop III stars are thought 
to have been more massive than typical stars forming today. 
The details of their formation process, especially the shape 
of the Pop III initial mass function (IMF), are still a sub- 
ject of intense study. Early works suggested these stars were 
extremely massive, exceeding 100 Mq, and formed one per 
~ Bromm et "aL"1|2002| [Yoshida et 



minihalo (Abel et al. 2000 



al. 2006 1. More recent studies, exploring fragmentation at 
higher densities ( |Turk et aL"1|2009| |Stacy et aL]|2010| |Clark 
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et al.||2011a"|b| |Greif et al.]|2011| |2012|, or modeling the ef- 



fects of protostellar feedback (Sta cy et al.|2012 Hosokawa et 
|al.|2011| ), have be gun to suggest an IMF extending to lower 
masses. Detailed simulations that will further constrain the 
Pop III IMF are needed to assess the role that these stars 
played in early cosmic milestones, such as chemical enrich- 
ment of the IGM, reionzation, and the formation of super- 
massive black hole seeds (e.g.JBarkana & Loeb 2001 Bromm 
fc Larson|[2004l |Ciardi fc Ferrara||2005| |Bromm fc Yoshida 
20111. Furthermore, determining the properties of Pop III 



stars is necessary for interpreting increasingly detailed ob- 
servations of high-redshift sources (e.g., |Dunlop||2012[) and 
of stellar relics in the local Universe (e.g., Frebel et al. 2005 



Beers fc Christlieb|2005l |Frebel fc Bromm|2010||Karlssonet1 



al.|2011| >. 

Directly observing a chemically pristine stellar popula- 
tion would represent a significant step towards understand- 
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ing the formation and properties of Pop III stars. How- 
ever, the chances of detecting individual, high-redshift Pop 
III stars with the upcoming James Webb Space Telescope 
(JWST) are very low (e.g. |Gardner et al.|2006| |Greif et al 



2009 Rydberg et al. 20101. While a single pair-instability 



supernova (PISN) from a massive Pop III star may be de- 
tectable (e.g.,|Wise fc Abel||2005| |Hummel et al.||2011| |Pan 



et al.||2012| |Tanaka et al.||2012[ ), clusters of Pop III stars, if 



they exist, would present the best opportunity for directly 
observing a metal- free stellar population ( Inoue|2011 Zack- 



risson et al.|2011[ ). It would seem that the haloes capable of 
hosting these clusters, with virial masses ~ 10 7 -10 8 M Q and 
often dubbed to be the first galaxies ( Bromm et al.|[2009 l, 
would have formed in high-density, biased regions, chem- 
ically pre-enriched with Pop III supernova ejecta ( Trenti| 



|& Stiavelli||2009 



Greif et al. 2010 Wise et al. 20121, pre- 



cluding the possibility of Pop III star formation. There are 
scenarios, however, in which these haloes reach the condi- 
tions necessary for atomic cooling while still metal-free. For 
example, if a strong hydrogen molecule dissociating 'Lyman- 
Werner' (LW) radiation background was set up sufficiently 
early (e.g., Haiman et al. 1997, M achacek et al.|2001 John 



son et aL]|2008 l or if Pop III stars ended their lives by col- 
lapsing directly to form black holes ( |Heger et al.|2003| ) , the 
onset of local metal enrichment would have been substan- 
tially delayed. Additionally, pockets of metal-free gas could 
have remained until very low redshifts owing to inhomoge- 
neous metal dispersal ( Scannapieco et al.|2002 Furlanetto fc| 
Loeb|2005||Tornatore et al.|2007||Trenti et al.|2009HStiaveHi| 
fc Trenti|2010| |Maio et al.||2010[ ) . Recently, the detectio n of 
metal-free gas clouds at z ~ 3 (Fumagalli et al. 20111 has 



confirmed that regions of space can remain chemically pris- 
tine long after reionization. 

It has been suggested that two different modes of metal- 
free star formation occurred in the early Universe. First gen- 
eration Pop III stars (Pop III. 1) formed from initial condi- 
tions completely unaffected by previous star formation. Sec- 
ond generation Pop III stars (Pop III. 2) formed from gas 
significantly influenced by the radiative output of previous 
star formation, but still containing no stellar nucleosynthetic 
products ( |0 'Shea et al.|2008||McKee fc Tan|20"08]|Bromm et 



al.|2009 1. Cooling by the hydrogen deuteride (HD) molecule 



is generally thought to differentiate Pop III. 1 and Pop III. 2 
star formation. Unlike H2, which cannot cool gas below 
~ 200 K, HD possesses an intrinsic electric dipole moment 
and can thus act as an effective cooling agent below 200 K, 
possibly resulting in stars with lower characteristic masses. 
The abundance of HD can be enhanced in regions with an el- 
evated free electron fraction. These regions can be produced 
from virialization shocks in haloes with virial temperatures 
T vir > 10 4 K (e.g.,|bh fc Haiman|2002[|Greif fc Bromm|2006 1 



or in the collapse of relic HII regions produced around Pop 
III.l stars dFerrara|1998l |Oh fc Haiman|2003||O'Shea et al 



2005 Yoshida et al.||2007[ ). For reference, T v j r is related to 
the virial mass of a halo as (e.g., Barkana fc Loeb|2001 ) 



2 x 10 



/ M viT 



V 10 8 m @ 



2/3 



1+Z 

10 



(1) 



where \i is the mean molecular weight {p, = 1.2 for neutral 
atomic primordial gas), mu is the mass of a hydrogen atom, 
fee is Boltzmann's constant, and M v i r is the total mass con- 
tained within the radius in which the average matter density 



is 187r 2 « 178 times the critical density. Theoretical investi- 
gations examining the chemistry, cooling, and dynamics of 
these regions have shown that gas is able to cool to the tem- 
perature of the cosmic microwave background (Tcmb) as a 
result of HD cooling. This may result in lower characteristic 
fragmentation masses (e.g., | Johnson fcT jromm 2 006[ ). 

One effect that can suppress gaseous collapse, star for- 
mation, and metal enrichment in small cosmic haloes is 
a pervasive UV background. Photons with energies in the 
range 11.2 eV < hp < 13.6 eV, the LW bands, are capable 
of photo-dissociating H2, the key cooling agent in metal- free 
gas below 10 4 K, and have a very small optical depth in a 
neutral IGM. Additionally, the photodestruction of H~, an 
intermediary in the gas-phase formation of H2, can also limit 
the H2 abundance. Many studies have explored the effect of 
an UV background on early structure formation ( Haiman ct 



al. 



al. 



1997 



2001 



Ciardi et al.|2000||Machacek et al.|2001||Ricotti et 
Mesinger et al.|2006||Wise fc Abel|2007| |Yoshidaet 



ah 2007 O'Shea fc Norman|2008| . It is accepted that above 
a certain radiation intensity, Jlw,2i ~ 10" 1 , LW radiation 
delays the collapse and cooling of metal-free gas until the as- 
sembly of more massive haloes Q Higher UV background in- 
tensities, Jlw,21 > 10, completely suppress baryonic collapse 
and cooling in haloes which allow only H2 cooling. In this 
regime, significant cooling will not occur until the assembly 
of larger mass haloes with virial temperatures T v i r > 10 4 K, 
the atomic cooling threshold. In these haloes, Lya emission 
will allow the gas to radiate its internal energy and col- 
lapse, effectively independent of the radiation background 
longward of 13.6 eV. 

Studies that have explored the thermodynamical evolu- 
tion of metal-free gas exposed to a strong UV background 
dOmukai||2001| |Omukai et al.||2008| |Safranek-Shrader et al.| 



2010||Wolcott-Green et al.|20il||Latif et al.|2011[ ) have sim- 
ilarly found that the evolution of metal-free gas undergo- 
ing free-fall or isobaric collapse is determined in large part 
by the background intensity, with the onset of effective H2 
cooling delayed with an increasing Jlw,2i- Gas that reaches 
the atomic cooling threshold is able to collapse isothermally 
while remaining at T ~ 8000 K until H2 forms in sufficient 
abundance for its cooling rate to exceed the adiabatic heat- 
ing rate, though this picture depends on the role played by 



Lya radiation trapping (see Latif et al. 20111. Addition- 
ally, there exists a spectrum-de 



Dendent critical LW radia- 



Omukai 



2001 



Shanget al. 20101, J£w,2i, 



tion intensity (e.g 
above which H2 never becomes an effective coolant. In this 
regime, LW irradiated gas collapses to high density via Lya 
and H~ free-bound emission. This evolutionary track has 
been suggested as a potential mechanism for the formation 
of supermassive black hole seeds via direct gaseous collapse 



(e.g.,|Bromm fc Loeb|[2003l I Begelman et al.|2006||Regan fc 
Haehnelt|2009| |Shang et al.|2010[).~ 



The atomic cooling threshold is an appealing criterion 
for classifying objects as the first galaxies (for a recent re- 
view see Bromm & Yoshida 2011 1. These haloes, with virial 



1 Here, Jlw,21 denotes the radiation intensity at the centre of 
the LW bands, 12.4 eV, in units of 10 -21 ergs -1 cm -2 Hz sr -1 . 
This is not to be confused with J21, the radiation intensity at the 
Lyman limit, 13.6 eV, in the same units. In general, Jlw,21 = 
J21 , where /3 = 3 fo r a 10 4 K blackbody spectrum an d 0.9 for a 
10 s K spectrum (e.g., Wolcott-Green fc Haiman|2011 1 
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masses > 5 x 10 7 Mq at z ~ 10, distinguish themselves Irom 
minihaloes in that in them, metal-lree gas can cool and col- 
lapse even in the presence of a strong UV radiation field. 
Additionally, supersonic turbulent gas flows, typically not 
present in minihaloes, should potentially develop during the 
assembly of the more massive haloes ( |Wise fc Abel|[2007| 



Greif et al."||2008| |Prieto et aT1|2012| >. This turbulence may 



have influenced the process of star formation in Pop III. 2 
star forming haloes. 

To understand star formation in the first galaxies, it is 
instructive to take guidance from the better understood case 
of star formation in the nearby Universe. Overall, the forma- 
tion of stars is observed to be extremely slow, in the sense 
that molecular clouds undergoing free-fall collapse should 
have star formation rates ~ 100 times higher than the ob- 



served rate (e.g., Zuckerman & Evans 1974 Krumholz & 
Tan|2007| |Evans et al.|2009||Kennicutt fc Evans|2012[ ). Ob- 
servational and theoretical work has suggested this ineffi- 
ciency stems from a number of effects, including supersonic 
turbulence, protostellar outflows, magnetic fields, and radia- 
tive feedback from massive stars. State-of-the-art radiative 
transfer simulations of star formation including these effects 
are able to reproduce the IMF and star formation rate in 



Orion- type Galactic star forming regions (e.g., Krumholz et 
aL]|2012| . While there are undoubtedly many effects which 
influence star formation, it is becoming accepted that the in- 
terplay of supersonic turbulence, gas self-gravity, and proto- 



stellar feedback are the key players (e.g., Mac Low & Klessen 



2004 Elmcgreen fc Scalo||2004| McKee fc Ostriker 20071, 



setting both the rate of star formation (e.g., Krumholz & 



McKee|2005l, and establishing the shape of the stellar IMF 


(e.g., Padoan & Nordlund |2002 


Padoan et al. 2007 Hen- 


|nebelle & Chabrier| |2008| |2011 


1. With this in mind, we 



are particularly interested in whether supersonic turbulence 
plays a role in regulating star formation in the first, metal- 
free galaxies. 

In this work, we present the results of a high-resolution 
cosmological simulation that follows the assembly of a metal- 
free, atomically cooling 3 x 10 7 Mq dark matter halo in 
an environment exposed to strong LW radiation. We ac- 
curately compute the column density of H2 to properly 
model the transition where H2 starts shielding itself from 
LW radiation, as is necessary for the formation of a cool, 
dense, baryonic core where star formation can take place. 
We utilize a spatially adaptive grid to resolve densities up 
to n = 10 s cm -3 and length scales down to ~ 1000 AU. 
Then, we employ sink particles to study the long-term frag- 
mentation tendencies of the gas. Similar studies that focused 
on LW suppression of H2 explored smaller values of the ra- 
diation intensity, Jlw,2i < 1, than we consider here and 
argued that H2 self-shielding was not an important effect 
in the small mass haloes they considered (e.g. 



Machacek 



et al.||200l] |Q'Shea fc Norman|[2008[ ). Other studies that 
explored much higher values of Jlw,21, relevant for the the- 
oretical scenario where supermassive black hole seeds form 
by direct gaseous collapse, only included H2 self-shielding in 
an approximate fashion based on purely local estimates of 



low J£w 21 an d thus permitting H2 cooling to play a major 
thermodynamic role. 

We will address two primary questions in this study. 
First, what is the nature of the fragmentation which occurs 
in the cold, self-shielding gas? Utilizing sink particles, we can 
evolve the gas long past the initial gravitational collapse to 
times when a significant mass has been accreted by these 
particles. We can also identify any additional collapsing re- 
gions which would represent other potential star forming 
clumps. In addressing this question, we will comment on the 
importance of HD cooling which is thought to give rise to en- 
hanced fragmentation, thus producing a distinct population 
of metal-free stars with lower characteristic stellar masses. 
Second, is it possible that LW radiation delayed collapse 
in metal-free atomic cooling haloes can produce clusters of 
Pop III stars that have luminosities high enough to be de- 
tectable and identifiable with the JWST? Zackrisson et all 



( 2011 1 showed that clusters of Pop III stars with total stellar 
masses as low as ~ 10 5 Mq should be detectable at z ~ 10 
in deep JWST exposures and that their primordial composi- 
tion can be ascertained based on simple colour criteria. In a 
similar study, |Inoue| ( |2011[ ) suggested that a star formation 
rate of a few solar masses per year at z ~ 10 will be needed 
for JWST detection. We provide rough estimates for the 
star formation efficiencies and mass spectra of high-redshift, 
metal-free stellar populations in atomic cooling haloes and 
comment on the feasibility of detection. 

The outline of the paper is as follows. In Section [2] we 
describe the initial conditions of the simulation and our nu- 
merical methodology. In Section [3] we present results of the 
simulation. In Section|4]we comment on the nature of super- 
sonic turbulence. In Section[5] we discuss the trends towards 
gravitational fragmentation of gas and attempt to predict 
the properties of the expected starburst. In Section[6]we dis- 
cuss the role of HD cooling. In Section [7] we provide further 
comments on our results, including a discussion on the ex- 
pected intensity of the LW background, the effect of internal 
radiative feedback, and the detection prospects of metal-free 
stellar clusters with the JWST. And finally, we summarize 
our results and conclude in Section [8] 

Throughout this paper we assume cosmological param- 
eters consistent with the Wilkinson Microwave Anisotropy 
Probe (WMAP) 7 year results ( TKomatsu et al.|2011[ ): Q A = 
0.725, n h = 0.0458, fi m = 0.275, h = 0.704, <r g = 0.810, and 
n s — 0.967. Additionally, all quantities will be expressed in 
physical rather than comoving units unless explicitly stated 
otherwise. 



2 NUMERICAL SETUP 

2.1 Algorithms and Initial Conditions 

We use the publicly available adaptive mesh refinement 
(AMR ) code FLASH ( |Fryxell et al.|2000||Dubey et al.|2008 



2009), version 3.3, which solves the equations of Eulerian hy- 
drodynamics with the directionally split, piecewise parabolic 



the H2 column density (e.g., Bromm & Loeb 2003 Shang method of Colella & Woodward (19841. Baryons are repre 



et al. 20101. This work bridges the gap between these two 
extremes, exploring the scenario in which collapse in metal- 
free gas is delayed by LW radiation until the halo reaches 
the atomic cooling threshold, but with Jlw,2i remaining be- 



sented by a multispecies fluid and dark matter by collision- 
less, massive particles. The gravitational potential of gas 
and dark matter is computed with the iterative multigrid 



Poisson solver of Ricker ( 2008 \ 
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We initialize the simulation at z — 146 in a 1 Mpc 3 (co- 
moving) box. Cosmological initial conditions were generated 
with MPGRAFIC ( |Prunet et al.|2008[ ), a parallel version of 
the multiscale Gau ssian random field generator GRAFIC 
( Bertschinger|p001 1 . We first run a 128 3 dark matter only 
simula tion and use the halo finder HOP of lEisenstein fc Hut I 
( 1998 1 to locate the site of the first 10 s Mq dark matter 



halo in the simulation volume. We then carry out a hier- 
archical zoom-in procedure to increase the mass resolution 
around this halo using three separate levels of dark matter 
refinement to reach a maximum effective resolution of 512 3 
and an effective dark matter particle mass of 230 Mq in the 
target halo itself. We choose the volume of the highest res- 
olution region such that the total mass contained within it 
is 1O 9 M0, 10 times the mass of our target halo. We have 
verified that only the highest resolution dark matter parti- 
cles are found in our target halo. Given our box size, the 
expected number of 10 8 Mq dark matter haloes at z = 10 is 
of order 10 (e.g., |Greif et al.|2008 |. 



2.2 Resolution and Adaptive Refinement Strategy 

In order to capture gaseous collapse to progressively higher 
densities, we utilize AMR which creates more finely spaced 
grids (refines) in localized regions. To trigger refinement, we 
use two separate criteria that are based on local gas prop- 
erties. We employ a criterion very similar to that used in 



Wise & Abel (20071 which refines based on baryonic over- 



density. In this scheme, the threshold comoving density for 
refinement is 

p th = 3p h 2 3 «- l <X 1+ V , (2) 

where pb = 3//of2b/87rG is the comoving baryonic density, 
{ is the current level of refinement, U is the initial level of 
refinement (given our grid resolution, k = 5), and 4> is the 
Lagrangian refinement factor; (j> = enforces constant bary- 
onic mass-per-cell while 4> < implies that the mass-per-cell 
decreases with increasing refinement level. We set <j> = —0.3 
which results in a baryonic mass-per-cell at the highest re- 
finement level (/ max = 22) of « 0.1 M©. 

A hydrodynamical simulation of a self-gravitating fluid 
must properly resolve the Jeans length to be physically re- 
liable. For grid-based Eulerian codes, this requirement is 
expressed by the Truelove criterion ( jTruelove et al.||l997[ ), 
which states that the Jeans length, 

/ 2\ 1/2 

Mi) (3) 

must be resolved by at least 4 grid cell widths to avoid arti- 
ficial fragmentation. While this criterion was originally for- 
mulated for isothermal gas and does not take into account 
the effect of Hubble expansion, it is commonly utilized in 
hydrodynamical-cosmological simulations. It would be un- 
necessary to enforce this criterion across the whole grid, 
thus we only apply it in the innermost region of the sim- 
ulation where the most highly refined dark matter particles 
are present. In this work, we always resolve the Jeans length 
by at least 12 grid cells and derefine the grid if it is resolved 
by more than 24. While this is more than sufficient to prop- 
erly resolve the fragmentation tendencies of gas, it may be 
insufficient to study the possible small-scale turbulent flow 
in the gas (see Federrath et al.|2011 1. 



The FLASH code allows for no explicit force soften- 
ing beyond local (e.g., one cell wide) cloud-in-cell smearing 
of the dark matter particle mass density on the computa- 
tional grid. In our simulation, because of AMR, the grid 
spacing can become many orders of magnitude smaller than 
interparticle separation. When this happens, the dark mat- 
ter particle discreteness can introduce severe artifacts into 
the calculation of the gravitational potential and cause gas 
on the computational grid to feel the gravity of individual 
dark matter particles. To achieve sufficient smoothness in 
the dark matter particle mass distribution, we have devel- 
oped an algorithm that spatially smears the dark matter 
density before it is passed to the Poisson solver. 

For each dark matter particle we compute a smoothing 
Kernel radius r s over which the mass of the particle is to be 
distributed, 



r s = 0.3 (A/ DM /p b ) 1/3 



(4) 



where Mum is the dark matter particle mass and pb is the 
baryonic density of the cell containing the dark matter par- 
ticle. This choice of smoothing radius ensures that the av- 
erage dark matter mass per cell contributed by a single 
dark matter particle inside its smoothing kernel, if multi- 
plied by Qb I (fim — Ob) , is not larger than the baryonic mass 
inside the particle's host cell. Within the smearing radius, 
we distribute the particle mass following the quadratic (or 
'Epanechnikov') kernel oc 1 — (r/r a ) 2 . 

In practice, we achieve the smearing by replacing the 
dark matter particle with a sufficient number of daughter 
particles with total mass equal to that of the parent particle 
and with density approximating the Kernel profile. Given 
the highly parallel nature of this simulation, this algorithm 
takes advantage of the FLASH code's capability for moving 
Lagrangian data between structured, Eulerian blocks (see 
Dubey et alj]2011| [2012| ). The daughter particle spacing is 
approximately equal to the cell spacing of the parent par- 
ticle's computational cell. Computation of the gravitational 
potential is performed from the daughter particle density. 
Computation of the parent particle's acceleration is then 
carried out in a momentum conserving fashion by summing 
up the gravitational force over daughter particles and as- 
signing the result to the parent particle. 

2.3 Chemistry 

Detailed, non-equilibrium chemistry is necessary to prop- 
erly model the thermodynamic state of cosmological gas 



flows (for a recent review see Glover 20111. Our chem- 
ical model tracks the most relevant chemical species in 
metal-free gas: H, H", H+, e", H 2 , Hj, He, He+, He ++ , 
D, D+, and HD. We evolve the species' abundances and 
the internal energy of the gas by simultaneously integrat- 
ing Aspects + 1 differential equations with a Bulirsch-Stoer- 



type, semi- implicit extrapolation mid-point method (Bader 



& Deuflhard 19831. We set our chemical timestep, which 



subcycles within the hydrodynamic timestep, as At = 0.1 x 
min{n e /|n e |,nH 2 /|nH 2 |,nHD/|nHD|}, where is the num- 
ber density of species i. 

The initial number density of hydrogen nuclei is 
given by n = pi,(z;)/[m H (1.0 + 4.0 zhc)], where p b {z) = 
(3#o/87rG) n h (l + z) 3 is the average, physical baryonic mat- 
ter density at redshift z, and ihc = 0.08 is the primor- 
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dial number fraction of helium. We set xu 2 = 2 x 10 , 
x H + — 3.8 x 10 -4 , and id = 4.3 x 10~ 5 , where Xi is the 
number density of species i relative to n, the number den- 
sity of hydrogen nuclei (henceforth the abundance of that 
species) . The abundance of electrons is calculated by enforc- 
ing charge neutrality. The initial gas temperature is set by 
assuming adiabatic cooling due to Hubble expansion after 
z ss 200, when gas and the CMB thermally decoupled. 

Properly modeling the formation of H2 and HD is es- 
sential since these molecules are the only low temperature 
(< 10 4 K) coolants in metal-free gas. In the absence of dust, 
H2 forms primarily through the gas-phase reaction mediated 
byH", 



H + e~ 
H~ + H 



H" + 7, 
H 2 + e~, 



(5) 



which is feasible in primordial gas owing to the residual ion- 
ization fraction present after recombination. At higher den- 
sities, > 10 s cm -3 , H2 can also form directly through 3-body 
reactions ( Palla et al.||1983 l, though since we do not simu- 
late densities this high, 3-body H2 formation will not be 
significant here. 

The HD molecule can be a significant coolant at tem- 
peratures < 200 K that can potentially cool the gas to Tomb- 
In the absence of a LW background, the HD abundance is 
primarily determined by its main formation pathway 



H 2 +D+^HD + H+, 



(6) 



where the D + to H + abundance ratio is set by the charge 
exchange reaction 

H + D + ^D + H + . (7) 

Given equilibrium in reacti ons |6| ) and Q, the HD abun- 
dance is (e.g., Omukai et al.|[2005[ ) 

xhd « 2exp(421K/T)a;H 2 x D . (8) 

Therefore with HD in chemical (not photodissociation) equi- 
librium, the HD to H2 abundance ratio can exceed the cos- 
mological D to H ratio by a large factor when T <C 421 K. 
This HD to H2 fractionation, which is due to the slightly 
higher binding energy of the HD molecule compared to that 
of H2 (Solomon & Woolf 19731, is a reason why HD is 



thought to be a significant coolant in Pop III. 2 star for- 
mation where additional cooling from H2 can induce an ele- 
vated HD abundance. As we shall show, even in the presence 
of a strong LW background, Equation (JsJ) describes the HD 
abundance fairly accurately, although since only a very small 
fraction of the gas cools down to T <C 400 K, significant HD 
to H2 fractionation does not occur. 



2.4 Gas Cooling 

Atomic and molecular radiative cooling processes have been 



well studied in astrophysical settings (e.g., Shapiro & Kang 



1987| |Cen||1992| |Osterbrock fc Ferland||2006[ ). For the den 

sities, temperatures, and chemical compositions relevant to 
this work, the most important cooling mechanisms are Lya, 
or more generally atomic, emission from neutral hydrogen, 
ro- vibrational emission from molecular hydrogen, and emis- 
sion from hydrogen deuteride. For completeness, we also in- 
clude Compton heating and cooling due to electron scat- 



tering of CMB photons, H and He recombination and colli- 
sional ionization cooling, and free-free emission. We do not 
consider cooling by metals or dust as we specifically focus 
on primordial gas. 

Lya cooling becomes significant above ~ 10 4 K when 
the electron fraction and gas temperature become high 
enough to excite this transition. This extremely efficient 
cooling channel can cool the gas to ~ 8000 K. Ro-vibrational 
emission from molecular hydrogen can potentially cool the 
gas further. Being a symmetric molecule, however, H2 lacks 
a permanent electric dipole moment and consequently can- 
not cool the gas below ~ 200 K. HD, if present, can cool 
gas to even lower temperatures, < 100 K, as it does possess 
an intrinsic dipole moment and more closely spaced energy 
levels. Nevertheless, an elevated electron fraction is still re- 
quired for it to form in significant quantities (e.g 



Johnson 



fc Bromm||2006[ ). We adopt the H2 cooling rate from Galli 
fc Palla| ( |1998[ ) and the HD cooling rate from |Flower et alT 
(20001). 



Finally, as the CMB imposes a lower limit on the tem- 
perature to which gas can radiatively cool, we adopt an ef- 
fective cooling rate of the form A c s(T) = A(T) — A(Tcmb), 
where Tomb = 2.7 K (1 + z) and A(T) is the total volumetric 
cooling rate not taking into account the CMB. This formu- 
lation ensures the gas temperature will not fall below Tomb 
unless it does so via adiabatic expansion. 

2.5 Sink Particles 

We utilize sink particles to follow the global evolution of 
the gas after the first baryon dominated region inside a halo 
undergoes runaway gravitational collapse. Without sink par- 
ticles, the first collapse would drive the computational mesh 
to arbitrarily high refinement, which would impose a pro- 
hibitively short Courant-Friedrichs-Lewy (CFL) timestep 
on the entire simulation. Sink particles allow us to self- 
consistently simulate the formation of multiple gravitation- 
ally bound collapsed structures over many free-fall times by 
accreting mass from the grid and putting an upper limit 
on the gas density. This computational method was origi- 



nally introduced in Bate et al. ( 1995 I and has been used 



extensively in numerical simulations, both in AMR (e.g., 
Krumholz et al.|[2004l |Wang et al.|[20T0l [Federrath et al.| 



2010[ |Padoan fc Nordlund||2011| |Girichidis et al.||2 011[) and 



smoothed particle hydrodynamics (SPH; e.g., Bromm et al 
2002l|Bate et al.||2003| |Stacy et al-pHO) |Greif et al.||2011 
Clark et al.|2011a[ ). 

Our sink particle implementation is identical to the 
method introduced in Federrath et al. (20101. Special care 



is taken to avoid the spurious creation of sink particles that 
may not represent gravitationally collapsing regions. To this 
end, we utilize additional checks for the creation of a sink 
particle in addition to requiring that the gas density in a 
cell be greater than some threshold pthrcsh- This includes 
enforcing that the local velocity divergence is negative and 
that the region of collapse is both gravitationally bound and 
Jeans unstable. 

Once formed, sink particles are capable of accreting 
mass directly from the computational grid. Cells with cen- 
tres within a constant distance, the accretion radius r acc , 
from a sink particle and with a density p > pthrcsh are ex- 
amined. For these cells, the mass increment AM — (p — 
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Pthresh)AV is calculated, where AV is the cell volume. Pro- 
vided that the cell velocity is radially directed towards the 
sink particle and that the mass increment AM is gravita- 
tionally bound to the sink particle and the gas within r acc , 
the increment is deducted from the grid and added to the 
sink particle. 

The gravitational interaction between sink particles and 
gas is computed by direct summation. To avoid extremely 
large accelerations, we employ cubic spline softening (e.g., 
Price & Monaghan 2007) which decreases the gravitational 
attraction of gas-sink and sink-sink interactions within r < 
r so ft, where r is the separation between a cell centre and a 
sink particle or between two sink particles. 

In the simulation here, we set the sink particle accretion 
and softening radii to be equal, r acc = r so f t = 0.01 pc ~ 
2000 AU, which is 2.5 times the grid spacing at the highest 
level of refinement. Additionally, we set the sink particle 



creation density threshold to pthresh 
corresponding to n = 10 s cm -3 . 



2.2 x 10" ib gcm" 



2.6 Transport of H2-Dissociating Radiation 

The occurrence of a metal-free atomic cooling halo requires 
previous suppression of Pop III. 1 star formation in mini- 
haloes that could have polluted the atomic cooling halo 
with metals. As previously discussed, this suppression could 
have been produced by a molecule-dissociating radiative 
background produced by neighboring star-forming galaxies. 
From the initial redshift, we impose a constant LW back- 
ground incident flux with intensity J21 = 100 (corresponding 
to Jlw,21 = 90) onto each of the six faces of the computa- 
tional box. This ignores periodicity of the domain, which 
is acceptable given that our target halo is located near the 
centre of the box. Due to absorption by neutral hydrogen in 
the IGM, we set J„ = shortward of the Lyman limit. We 
assume the spectral shape of this source is a 10 K black- 
body, representative of massive, metal-free stars (Bromm 
et al. 2001 Schaerer 2002 1 . While it has been suggested 
that the photodissociation of H~ can be crucial in regu- 
lating the H2 abundance, we neglect it as it has been shown 
to have minimal importance for 10 5 K blackbody sources 
(e.g., Omukai et al. 2008). However, if Pop III stars were 
not extremely massive (~ 100 M@) but had more moderate 
masses (~ 10 Mq), the appropriate radiation source tem- 
perature would be somewhat lower. We note, though, that 
the primary consequence of including H~ photodissociation 
in this context would be a shift in the density at which H2 
cooling becomes effective, and would not significantly alter 
our overall results. 

The photodissociation of H2 occurs through Solomon 



process ( Stecher & Williams 1967 1 in which a photon with an 



energy coinciding with either the Lyman or the Werner (LW) 
bands of H2 places the molecule in an excited electronic 
state. The subsequent radiative decay has a ~ 15% chance of 
reaching the ground state continuum, which results in molec- 
ular dissociation. Accurately modeling this process requires 
detailed modeling of hundreds of LW lines which is compu- 
tationally unfeasible in 3D hydrodynamic simulations. For- 
tunately, the photodissociation rate of molecular hydrogen 
can be expressed as (e.g., Abel et al.|1997 l 




10" 1 10° 



shield, colutr 



Figure 1. Comparison of the H2 self-shielding factor, Equation 
{11} , between approaches where the H2 column density ATjj, is 
calculated via our column density method (/shield column) and 
two common approximations (/shield, approx)- In these approxi- 
mations, Afjj 2 is calculated locally with a Sobolev (green squares) 
and a Jeans length approach (blue diamonds). The dashed line 
represents a one-to-one mapping. The Sobolev approximation for 
Nn 2 produces an H2 self-shielding factor in close agreement with 
our method, except for a few cases of disagreement, likely due to 
small velocity gradients. A hybrid of these two approaches may 
be useful in future work (see |CIark et al.|2011afl . 



where the dimensionless factor /shield, h 2 ^ 1 accounts for 
H2 self-shielding with /shield, h 2 = 1 corresponding to no 
shielding and /shield, h 2 = to complete shielding. For a 
static, cold medium, /shield. h 2 can be written solely as a 



function of the H2 column density A r H 2 (Draine & Bertoldi 
fT996l, 



shield, Ho = mm 



1.0, 



H 2 



10 14 cm" 2 



(10) 



To more accurately model a dynamic medium, |Draine fc| 
|Bertoldi| (|1996[) also provided a fit for /shield, h 2 which takes 



(11) 



thermal gas motion into account, 
0.965 



/shi 



cld,H 2 = 



(l + x/b 5 )- 
0.035 

+ (l + a;)0.5 



x exp[-8.5 x 10" 4 (1 + . 



a = 2, 



where x = iVH 2 /5 x 10 14 cm~ 2 , b 5 = b/10 5 cms~ 
and b = 9.12kms~ 1 (T/10 4 K) 1/2 is the v elocity spread pa- 
rameter for H2 (e.g., Ahn & Shapiro 20071. Unless otherwise 
noted, we will use Equation (111 with a = 1.1, a modifica- 
tion suggested by | Wolcott-Green et al.| ( |2011[ ), to calculate 
the self-shielding factor for H2 from LW radiation. 

Hydrogen deuteride can also be destroyed by LW radi- 
ation at a rate similar to that given in Equation {9|, 



fc HD = 1.5 x 10 _i V L 



- /shie 



I /shi 



cld,H 2 ,HD S 



(12) 



which includes both a factor due to HD self-shielding and 
a factor accounting for the shielding of HD by H2. The HD 



self-shielding factor is equivalent to that in Equation (111 



kn 2 



1.38 X 10 JlW, 21 /shield, H 2 S 



(9) 



with the HD column density Nud replacing Nu 2 ■ Due to 
slight energy differences between the H2 and HD LW line 
centres, the shielding of HD by H2 does not effectively occur 
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until there is a large H2 column density, 7Vh 2 



Wolcott- Green & Haiman (20111 provided a fit for the HD 



x 10~ 3 x) 



(13) 



shielding factor due to H 2 

/shicld,H 2 ,HD = ^ + ^o. 2 3 8 ex P (-5.2 

where x = Nn 2 /2.34 x 10 19 cm" 2 . Both H 2 and HD can also 
be shielded from LW radiation by neutral hydrogen, however 
we neglect this effect as self-shielding is always the dominant 
effect. In fact, HD photodissociation is never significant in 
determining the HD abundance. Instead, as we will argue, it 
is primarily the H2 photodissociation and self-shielding that 
determine both the H2 and HD abundances. 

We compute the molecular column densities, Nh 2 an d 
JVhd, using an on-the-fly, non-local approach very similar 
to the 'six-ray' approxim ation (e.g., |Nelson fc Langer|1997| 



Glover fc Mac Low|2007| |Glover et al.|2010[ |. Specifically we 
compute the column density at each point in our box as 



Nn 2 (r)=min j {N H2 j(r)} 



(14) 



where j — ±x, ±y, ±z, and Nn 2 ,j is the column density from 
point r to the edge of the box along direction j. The same 
approach is used for calculating the HD column density. The 
presence of bulk velocity gradients, which our approach does 
not account for, may act to decrease the importance of self- 
shielding by Doppler shifting absorption line centres (e.g., 
Draine fc Bertoldi|[l996} |Glover fc Brand|[200T| ), although 
Wolcott-Green et al. 1 2011 1 argue that this effect is minimal 



in radially coherent gas flows. 

While we self-consistently compute the column densi- 
ties to each point in our simulation, a commonly used (e.g., 



Bromm & Loeb 2003 Shang et al. 20101, and less com- 
putationally expensive, approximation is Nn 2 « nn 2 L c h a r, 
where Lchar is a local characteristic length scale and nn 2 
is the number density of molecular hydrogen. In Figure [I] 
we compare /shield, h 2 computed using Nn 2 calculated with 
our six-ray approach to / s hi c id,H 2 computed with Nn 2 based 
on two local characteristic lengths L c har at the time when 
the gas density first reaches 10 8 cm -3 . The first is the Jeans 
length, I/char = Lj, which is motivated by the assumption 
that the majority of shielding occurs in a region of size sim- 
ilar to the local Jeans length. The second approximation is 
a close analogue to the Sobolev length ( Yoshida et al.|2006 
Clark et al. 2011a I, L c har = c s /|V ■ v\, which, assuming a 



constant velocity gradient, is the distance at which absorp- 
tion line centres are Doppler-shifted by approximately one 
thermal line width. As is evident in Figure [I] both approxi- 
mations agree reasonably well with the more accurate, non- 
local approach, with the Sobolev approach yielding slightly 
better overall agreement. While this agreement is likely co- 
incidental given the disparate physics involved in all three 
approaches, it is still reassuring that they are all of similar 
magnitude. At least in this regime, a combination of the two 
local approaches could be useful in future work (e.g., [CTark| 
et al.|2011a\ 



3 RESULTS 

The LW radiation field of intensity J 2 i = 100 prevents the 
H2 abundance from reaching the level that would permit 
efficient gas cooling in any halo not capable of atomic line 




10° 10 1 
Radius [pc] 



icr 



Figure 2. Radial profiles of average gas density centered around 
the point of maximum density with time measured from the 
point of sink particle formation which occurs shortly after the 
gas first reaches n = 10 8 cm -3 . The gas density always maintains 
a roughly p oc r~ 2 density profile (shown by the straight line) with 
departures from the formation of the core and the gravitational 
influence of the sink particle at the final time shown. 

cooling]^] Hydrostatic collapse still occurs in these haloes, 
but only along an adiabat. At a redshift z ~ 13, our target 
halo can be classified as an atomic cooling halo, with an 
average gas temperature of 8000 K and most of gas within 
the virial radius lying between 2000 K < T < 1.2 x 10 4 K. 
With gas now able to radiate its internal energy, collapse 
proceeds isothermally at T ~ 8000 K resulting in an n oc r~ 2 
density profile within the virial radius — see Figure [2] 

When the maximum gas density near the halo centre 
reaches ~ 200 cm" 3 at z = 12.1, the virial mass of the sur- 
rounding halo is 3 x 10 r M© with a virial radius « 750 pc, 
just nominally fulfilling the standard analytic atomic cool- 
ing criterion, Equation |T]). At this point, the H 2 cooling 
rate first exceeds the rate of adiabatic heating, resulting in 
rapid gas cooling to T ~ 400 K. It should be noted that 
H2 self-shielding is not the determining factor in setting the 
density at which H 2 cooling becomes effective. Once that 
happens, however, the self-shielding factor rapidly drops 
below unity indicating strong shielding. Indeed, the inclu- 
sion of self-shielding is essential for the cooling instability to 
continue and for the H2 abundan ce to eventually reach the 
asymptotic abundance of ~ 10 -3 ( Oh fc Haiman|2002 1. 

In a suite of cosmological simulation exploring the ef- 
fect of LW radiation on the fraction of cool gas in haloes, 



Machacek et al. (2001 1 determined the threshold halo virial 



mass for metal-free gaseous collapse as a function of radia- 
tion intensity to be 



A/th = {1.25 x 10 5 + 2.9 x 



10 6 [J LW ,2i] a47 }M e 



(15) 



where we have converted their -Flw into Jlw,2i- Extrapo- 
lating beyond the intensity range, Jlw,2i < 0.080, explored 



by Machacek et al. (20011, Equation (151 still yields an ac- 
curate prediction for the halo mass M«3x 10 7 M© where 



1 If we included no radiation background, the first minihalo in 
our box would have collapsed at z ~ 19; see |Ritter et a,l .] ( | 20 1 2 1 ) , 
who used identical cosmological initial conditions. 
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Figure 3. Density-temperature phase plot at the time of sink 
particle formation at z = 12.1 showing gas within 2ii v i r of the 
halo. The colour corresponds to the mass contained in the phase- 
space region. The solid black curve is a one-zone calculation of 
a thermodynamic evolution of metal-free gas undergoing free-fall 
collapse with no radiation background (Pop III.l track). After 
the gas in the simulation is able to cool via H2 cooling starting 
around n ~ 2 X 10 2 cm~ 3 , most of the cooled gas reaches min- 
imum temperature at T ~ 400 K and n ~ 10 4 cm -3 . The gas 
that continues to collapse converges towards the standard Pop 
III.l molecular track. Adiabatic gas cooling due to expansion in a 
turbulent medium is apparent around n ~ 10 3 cm' 3 , permitting 
some gas to reach temperatures T < 100 K. 



self-shielding and collapse first occur in the simulation here. 



Since Machacek et al. (20011 did not run simulations with 



stronger radiation fields, their simulations did not produce 
any haloes in which atomic line cooling was effective. In- 
deed, once an atomic cooling halo forms, Lya cooling per- 
mits rapid isothermal collapse provided that the halo mass 
exceeds a minimum that depends very weakly on the LW 
radiation intensity. 

We will denote the density at which the H2 cooling rate 
becomes larger than the rate of adiabatic heating and the gas 
evolution leaves the atomic isothermal track with n coo i. We 
show a representative density-temperature phase diagram in 
Figure [3] which shows this density is n coo i ~ 200 cm" 3 and 
this density is found within a radius of ~ 10 pc. This is a 
larger density than n coo i « 10 cm -3 found by Shang et al. 
(20101 in a similar simulation with an identical radiation 



background. We attribute this discrepancy to our different 
treatments of H2 self-shielding, in particular the method for 



computing Nn 2 . Shang et al. ( 2010 1 utilized a local approx- 
imation A r H 2 = «h 2 Lj, while here we computed Nn 2 more 
accurately as described in Section |2.6| To further under- 



stand the source of this discrepancy, we can equate the H2 
gas cooling time i CO oi,H 2 and the free-fall time tg to obtain 
a rough estimate of n coo i, 



^cool 



180 cm" 



•i-n 2 



3 x 10- 7 



(16) 



where xn 2 is the H2 abundance at n coo i, iff = [37r/(32Gp)] 1/ ' 2 
is the free-fall time, i coo i,H 2 = 3nfcBr/2An 2 is the H2 cool- 
ing time, and Ah 2 is the volumetric H2 cooling rate (e.g., 
Galli fc Palla||1998[ ) which we evaluate at 8000 K. Equation 



( 16 1 is simply n coo i in terms of the required H2 abundance 
for the molecular cooling rate to exceed the rate of adia- 
batic heating, where 3 x 10~ r is the H2 abundance at which 
the transition to efficient molecular cooling takes place in 
our simulation. We can simplify this further by reasonably 
assuming that the H2 abundance is determined by its pho- 
toequilibrium value, 

kt u— x c 



*^H 2 ; photo — 



kv 



(17) 



where k { H - is the formation rate of H~, ku 2 is the pho- 
todissociation rate of H2 given by Equation |9|, and HP 
photodestruction is considered negligible. Inserting Equa- 
tion (171 into Equation (16 1 and assuming /shield, h 2 = 1 



results in 



f-cooi ~ 210 cm 



f J21 

Vioo 



2/3 



Xe 

5 x 10" 5 



-2/3 



(18) 



where J21 is the unattenuated radiation intensity and x e 
is the free electron abundance at n coo i. Equation (18 1 is in 



excellent agreement with the density at which H2 cooling 
becomes effective in the simulation, suggesting that self- 
shielding is not necessary for the onset of H2 cooling. Al- 
lowing for self-shielding, /shield, h 2 ^ 1, and adopting the 



prescription of | Shang et al. (20101 with Equation (10 1 and 
A r H 2 = wh, Lj gives us 



^cool 



10 cm 



-3 / J21 
\100 



-2/3 



(19) 



where we have now normalized to a slightly higher value 
of the electron abundance found in the simulation at lower 



densities. Equations ( 18 1 and ( 19 1 separately agree reason- 



ably well with our simulation and with that of |Shang et] 
al. (20101, respectively, suggesting that the H2 column den- 



sity computed with the local Jeans length is an overestimate 
compared to that computed with the more detailed radia- 
tive transfer described in Section [2.6| As we shall argue in 
Section [7. 3| this may have an impact on the star formation 
efficiencies and observability of these systems. 

The onset of effective H2 cooling at n coo i leads to the 
rapid emergence of a cold < 10 3 K, dense n ~ 10 4 cm -3 
core, similar to that forming in the process of Pop HIT star 
formation (e.g., [Bromm et al.||2002 l. While we have argued 
that self-shielding is not important in determining n CO oi, the 
H2 self-shielding factor does begin to drop below unity as 
the gas cools, reaching as low as /shield, h 2 w 10" 6 at the 
highest densities. Proper treatment of self-shielding is there- 
fore essential for accurately computing the chemical state of 
high density gas. To discriminate cold, molecule-rich, self- 
shielding gas from the warmer outer halo, we will generally 
utilize a simple criterion where we select the cells that have 
an H2 shielding factor /shield, h 2 < 10 -2 . Some of our results 
may be sensitive to how we decide which gas belongs to the 
self-shielding core, although we have verified that alterna- 
tive criteria, such as only selecting gas with temperatures 
< 2 x 10 3 K or densities > 10 3 cm~ 3 , would yield similar 
conclusions. 

After effective cooling begins when the gas density 
reaches n coo i, gas collapses at close to the free-fall rate and 
reaches n = 10 8 cm -3 in ~ 3 Myr. As the collapse progresses, 
we can estimate the mass of the first gravitationally unstable 
clump by comparing the enclosed gas mass to the Bonnor- 
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Figure 5. Mass-weighted line-of-sight gas density projections 3.5 X 10 5 years after the sink particle formation (z = 12.1). The six panels 
show progressively smaller fields of view. The black dot marks the location of the sink particle while the black circle in the top panels is 
the virial radius (~ 750 pc) of the target halo. The spatial scale, in physical units, is shown at the bottom of each panel. The upper-left 
panel shows neighboring haloes and the clustered cosmological environment where our target halo formed. Self-shielding, cold gas is seen 
in the bottom-left panel and the filamentary, irregular nature of this gas is apparent. The high density disc which forms around the sink 
particle is visible in the bottom-right panel, approximately face-on. 
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Figure 4. Accretion rate (top panel) and mass (bottom panel) 
of the sink particle as a function of time. Power law fits to the 
accretion rate and mass in the first 10 5 years are shown in red. 
The time-averaged accretion rate (blue dashed line) of the sink 
particle is 3x 10 _3 Mq yr — , which we take to be an upper limit on 
the star formation rate that occurred in this time period. The late 
time accretion rate approaches 10 — 3 MQyr — 1 . We also show the 
characteristic accretion rate c 3 /G (green dot-dashed line) with 
the adiabatic sound speed c s evaluated at T = 800 K. 



Ebert mass (lEbert 1955 Bonnor 1956 \ at different radii 



The Bonnor-Ebert mass can be written as 



Mb 



pl/2 G 3/2 



(20) 



where ay = (k^T / ^.m-n) 1 / 2 is the isothermal sound speed, 
Po is the ambient pressure, and mi = 1.18 is the maximum 
dimensionless mass in the solution to the Lane-Emden equa- 
tion (e.g., |Stahler fc Palla|[2005l ). In Figure [6] we show the 
enclosed gas mass (solid line) and the Bonnor-Ebert mass 
(dashed line) both as a function of distance from the point of 
maximum gas density, approximately 10 s years before sink 
particle formation. To evaluate the Bonnor-Ebert mass as a 
function of radius, Mbe(t), we take Po to be the pressure at 
radius r and a-r to be the mass-weighted isothermal sound 
speed interior to r. The enclosed gas mass first exceeds the 
Bonnor-Ebert mass at a radius of ~ 1 pc, corresponding to a 
mass of ~ 10 3 M . We thus expect the first gravitationally 
unstable fragment to have roughly this mass, though there 
remains the possibility that additional fragmentation may 
occur on unresolved scales or at later times. 

Sink particles are allowed to form at a density of 
n = 10 s cm -3 to follow the gas evolution for multiple 
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Figure 6. Enclosed gas mass (solid line) and Bonnor-Ebert mass 
(dashed line) as a function of distance from the point of max- 
imum density roughly 10 5 years before sink particle formation. 
The maximum radius where the enclosed gas mass exceeds the 
Bonnor-Ebert mass provides a rough measure of the mass of the 
first gravitationally unstable clump. We see the onset of gravita- 
tionally instability occurring inside a radius of ~ 1 pc when the 
enclosed gas mass is 1O 3 M . 



free-fall times (see Section 2.5 I. We run the simulation for 



3.5 x 10 yrs after the formation of the sink particle. At 
this point, the self-shielding core has an average density 
r;8x 10 3 cm -3 , total mass w 1.7 x 10 4 Mq, mass- weighted 
root-mean-square (rms) velocity n rms ~ 7.1 kms -1 , aver- 
age temperature m 480 K, and rms Mach number A4 = 
Wrms/Cs ~ 3.3. Additionally, the sink particle has increased 
in mass to « 1100 A/0, which is ~ 6% of the total mass 
of the self-shielding gas. No additional sink particles form, 
meaning that additional gravitational fragmentation, if any, 
would have been confined to occur on unresolved scales in- 
side the sink particle. 

In Figure |4j we plot the accretion rate and mass of the 
sink particle as a function of time. Initially, the sink par- 
ticle accretion rate is M s i n k ~ 0.1M©yr -1 which drops to 
0.01 M0yr -1 within 10 4 yrs. The time averaged accretion 
rate over the whole simulation is 0.003 Mq yr _ . Given the 
average density and mass of self-shielding gas at the end of 
the simulation, n « 10 4 cm -3 and « 2 x 1O 4 M0, respec- 
tively, the average rate of accretion onto the sink particle is 
roughly a factor of 10 lower than what would be expected 
if the gas had been undergoing free-fall collapse. We discuss 
this inefficiency of gaseous collapse further in Section [5] 

In Figure[5]we show mass- weighted line-of-sight gas den- 
sity projections on six different spatial scales 3.5 x 10 5 yrs af- 
ter the formation of the sink particle. The upper-right panel 
shows neighboring haloes and demonstrates the clustered 
cosmological environment where the target halo formed. The 
bottom-left panel displays the extent of self-shielding gas 
and its irregular, filamentary density structure. The self- 
shielding core is small, ~ 5 — 10 pc, compared to the virial 
radius of the halo, P v i r ss 750 pc. The bottom-right panel 
shows that a rotationally supported disc, ~ 0.1 pc in diame- 
ter, forms around the sink particle approximately 10 5 years 
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Figure 7. Abundance of H2 (top), HD (middle), and free elec- 
trons (bottom) as a function of density at the time when the 
sink particle formed. In the top panel, the dash-dotted line is the 
photo-equilibrium H2 abundance not considering self-shielding 
(/shield, h, = 1) while the solid line includes self-shielding. As 
is shown, self-shielding is essential for H2 abundances larger than 
~ 10~ 6 . In the middle panel, the black line is the equilibrium 
abundance of HD only assuming reactions involving H2 (see Equa- 
tion [sj and excluding HD photodissociation. The three equilib- 
rium curves were computed with data extracted directly from the 
simulation and represent mass-weighted averages as a function of 
density. 



after the sink's formation. The disc is seen approximately 
face-on. 

In Figure [7] we show the abundances of H2, HD, and 
free electrons as a function of density at the time of sink 
particle formation. In gas with a density of n — 10 4 cm -3 , 
the H2 abundance reaches ~ few x 10 -4 while the HD abun- 
dance is ~ few x 10~ 8 . Both of these values approximately 
match equilibrium abundances, assuming that H2 is in di- 
rect photo-dissociation equilibrium with the self-shielding- 
attenuated LW radiation field and HD is close to equilib- 
rium with H2 via the reaction in Equation Jfjl. We show 
the equilibrium abundance expectations for H2 and HD with 
solid black lines. Above this density, H2 approaches the non- 



equilibrium asymptotic abundance of xu 2 ~ 10 ' (e.g., Oh 



fe Haiman||2002 1 while HD remains in approximate equilib- 
rium with H2. We observe only minor fractionation of HD 
over H2, which is to be expected given the relatively high 
temperature of the gas ~ 400 K compared to what is needed 
for the xhd/xh 2 ratio to exceed the cosmological id/% ra- 
tio. In fact, the HD abundance is generally just below its 
equilibrium value at densities 10 cm~ 3 < n < 10 s cm . 
The electron abundance shows the expected behavior - 
when efficient H2 cooling begins at n ~ 10 2 cm -3 , the abun- 
dance is below the residual recombination value even though 



there is significant collisional ionization, 



10 , at lower 



densities (not shown in Figure[7|. We propose this as one ex- 
planation for the lack of significant HD cooling and discuss 
this further in Section [6] 



3.1 Kinematical State and Evolution 

To gain further insight into the properties of the gas flow we 
proceed to analyze the gas kinematics at different stages of 
the collapse. Before the onset of effective H2 cooling, Lya 
cooling acts as a thermostat that keeps gas inside the halo 
at a temperature ~ 8000 K. Towards the centre of the halo, 
isothermal collapse proceeds under quasi-hydrostatic con- 
ditions. However, with increasing density, the cooling time 
decreases faster than the free-fall time, eventually reaching a 
point at which H2 cooling is so rapid as to effectively remove 
pressure support and set the gas into free-fall collapse. This 
also marks the onset of self-shielding which allows the H2 
abundance and cooling-rate to increase further. This ther- 
mal instability proceeds until the gas temperature is roughly 
400 K. The collapse is also responsible for exciting some bulk 
turbulent motions and increasing the Mach number of the 
flow. At the time of sink particle formation, the self-shielding 
gas has an rms velocity « rms ~ 7kms _1 corresponding to a 
Mach number Ml ~ 3. 

In Figure [8] we show different components of the gas 
velocity computed in annular shells at four times: (a) just 
before the onset of effective H2 cooling, (b) after the onset 
of H2 cooling but before the maximum gas density reached 
10 8 cm -3 and a sink particle formed, (c) at the moment of 
formation of the sink particle, and (d) at the end of our sim- 
ulation, 3.5 x 10 5 years after sink formation. We consider five 
components of velocity, each computed via mass-weighting 
in annular shells centered on the sink particle location, ex- 
cept in panel (a) where it is centered on the point of maxi- 
mum gas density: 



(i) The rms velocity 
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Figure 8. Mass- weighted average gas velocities computed in annular shells. Panels (a) through (d) correspond to times —1.2 X 10 6 , —9.9 X 
10 4 , 7, and + 3.5 X 10 s yrs from the formation of the sink particle. We show the rms velocity minus the centre of mass motion of the halo 
(black solid line), sound speed (red triple-dot dashed line), radial velocity (orange dashed line), rms turbulent velocity (green dotted 
line), and rms rotational velocity (blue dot-dashed line). See text for additional details. 



'V^ I |2\ 1/2 



(21) 



where the sum ranges over all cells with centres inside the 
annular shell, m is the cell mass, and v is the gas velocity 
minus the centre-of-mass motion of the particular annular 
shell. 

(ii) The adiabatic sound speed c s = ('ykBT/ nmu) 1 ^ 2 , 
where T is the mass-weighted temperature inside the shell. 

(iii) The rms rotational velocity of the mean rotation in- 
side the shell 



Vrot 



I |2\ 1/2 



(22) 



where v Iot ,i = Vi x fi is the azimuthal rotational velocity of 
the mean flow in cell i, ft = L/I is the mean angular veloc- 
ity inside the shell, L = Y2i m i r i x v i is the total angular 
momentum inside the shell, and I = Ei m 'l r i x L\ 2 is the 
moment of inertia of the shell in the direction defined by the 
total angular momentum. 

(iv) The radial velocity of the bulk motion inside the shell 



«rad 



J2i TTlifi ■ Vi 



(23) 



(v) The rms turbulent velocity relative to the mean flow, 
quantifying unordered gas motions 



fturb = 



where v T 



J2i m i\ V i 



1/2 



frad Ti is the radial velocity of cell i. 



(24) 



As shown in Figure [8] before the onset of significant H2 
cooling (panel a), typical rms gas velocities in the halo are 
~ lOkms -1 near the virial radius and drop to ~ 5kms _1 
in the inner ~ 10 pc. The sound speed is greater than the 
rms velocity inside 300 pc, somewhat smaller than half of 
the virial radius. Inside the virial radius bulk gas motions 
are dominated by turbulent and rotational motions, as ra- 
dial infall motions greatly decrease after the virial shock. 
Gas kinematics inside the halo can be described as a tur- 
bulent, subsonic flow. Within 1 Myr after the onset of H2 
cooling (panel b), and shortly before sink particle forma- 
tion, the mass of self-shielding gas is w 10 4 Mq. The tem- 
perature has dropped to ~ 400 K in the inner ~ 1 pc. The 
rms velocities have increased within r w 10 pc where the gas 
flow is transonic, M ~ 1. At the time of sink particle for- 
mation (panel c), turbulent motions still dominate within 
~ 10 pc, although towards smaller radii the contributions 
of rotational and radial motions increase. The Mach num- 
ber is on the order of 2. Panel (d) shows the kinematical 
state of the gas 3.5 X 10 5 yrs later at the end of the simu- 
lation. Turbulent motion is still the dominant form of gas 
motion, particularly in the range 0.1 pc < r < 10 pc. Within 
r ~ 0.1 pc there is kinematical evidence of a disc with or- 
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dered rotational velocities approaching ~ 15kms _1 . The 
Mach number of the flow has also increased within the self- 
shielding core to M ~ 2 — 4. We will proceed to analyze the 
impact of these supersonic gas motions on the collapse of 
gas in the following section. 



4 SUPERSONIC TURBULENCE 

Supersonic turbulence is believed to play a key role in 
present-day star formation as it compresses gas to densities 
at which it can be susceptible to gravitational fragmenta- 
tion and collapse. It additionally can act as a source of pres- 
sure support against gravitational collapse on larger scales 
(e.g., Mac Low & Klessen 2004). Minihaloes, the sites of Pop 
III. 1 star formation, are not thought to contain fully devel- 
oped supersonic turbulence (e.g., Yoshida et al.||2006 l. The 
virialization process, gravitational inflow, and larger virial 
velocities of atomic cooling haloes are believed to give rise 
to supersonically turbulent flows, likely influencing star for- 
mation ( |Wise fc Aberj|2007| |Greif et al-1|2008| |Prieto et al~ 
|2011[). The turbulence which develops in the self-shielding 



core in our simulation is also generated by gravitational in- 
fall (see Federrath et al.||2011 for a study of general prop- 
erties of gravity-driven turbulence), but on a smaller spatial 
scale defined by thermal instability (e.g 
2002). 



Kritsuk & Norman 



The classical theory of Kolmogorov (Kolmogorov 1941 
Frisch 1995 1 describes turbulence as a chaotic fluid motion in 



which energy progressively cascades to smaller length scales 
through a series of eddies. These eddies eventually reach a 
small enough scale where viscous forces effectively dissipate 
the turbulence. Typical gas in the interstellar medium is 
highly compressible and compressible turbulence is known 
to result in a self-similar network of interacting shocks, large 



density contrasts, and a filamentary morphology (e.g., Krit- 
|suk et al.|2007| |. 

Many simulations have shown that driven, isothermal, 
supersonic turbulence results in a gas density probability 
distribution function (PDF) that is accurately described 
by a log-normal distribution (e.g., Vazquez-Semadeni|| 1994 



Passot fc Vazquez-Semadeni||1998| |Scalo et al.||1998| |Os 
triker et al.||2001 l. This log-normal shape is understood to 
be the result of multiple, independent shocks altering the 
logarithmic gas density contrast in random walk fashion, 
thus driving the density PDF towards a log-normal form. 
The addition of self-gravity, magnetic fields, or a realistic 
equation of state can alter this general shape. We write the 
lognormal distribution as 



p(s) ds - 



(27TO- 2 ) 



1/2 



exp 



ds , 



(25) 



where p(s) is the probability that a parcel of gas has s in 
the range [s, s + ds], s = ln(n/no) is the logarithmic density 
contrast, a is the standard deviation, and no is the average 
density. Additionally, by considering the normalization of 
the log-normal distribution, it can be shown that s — ~a 2 /2. 
Unless otherwise noted, we consider no/time-weighted quan- 
tities when discussing density PDFs. 

In Figure [9] we show the volume- weighted density PDF 
of self-shielding gas, defined such that /shield, h 2 < 10~ 2 , at 
five different times. As expected, the peaks of the density 




1 2 3 

l°9,o ( n /n ) 

Figure 9. Volume-weighted density PDFs of self-shielding gas 
(/shield, H 2 < 10~ 2 ) at various times from the creation of the 
sink particle. We show a log-normal fit to the final density PDF 
(black curve) which has a standard deviation a = 0.93. The PDF 
populates higher density contrast as gas collapses and at high 
densities shows a marked departure from the initial log-normal 
like distribution which can be attributed to the self-gravity of the 
gas. The dashed black line indicates a power law of slope —1.5, 
the expected PDF slope of an p oc r~ 2 density profile. 



PDFs lie just below the average density as most of the vol- 
ume is underdense, typical of compressive turbulence (the 
peak of a mass-weighted density PDF would lie just above 
the average density) . The black line is a log-normal fit to the 
final density PDF at t = 1.2 x 10 5 yrs, although such a fit 
to the PDF at any time would look approximately the same 
in this view. For this particular log-normal fit, a = 0.93. 

As time progresses, the gas density PDF develops 
a power-law tail towards higher density contrasts. These 
power-law tails are understood as a consequence of the self- 
gravity of the gas and have been seen in both numerical 
simulations of turbulence (e.g., |Scalo et 



al 



1998 



Klessen 



2000| |Federrath et al.|[2008| |Collins et al.||20il| |Kritsuk "el 
al.||2011| |Cho fc Kim||20l"T l and observations of active star 



forming complexes (e.g., Kainulainen et ah]|2009 Schneider 



|et al.|p012| . It can be shown that the slope of the den- 
sity PDF for a spherical p oc r~ n density profile is — 3/n 
(e.g., |Kritsuk et al-pHl"] |Federrath et al.||2011[ ). Thus, for 



an isothermal, n — 2, profile one would expect a power-law 
PDF with a slope of —1.5; this slope is shown in Figure [9] 
as a dashed black line. The slope of the power-law tail just 
after sink particle formation matches this expectation quite 
well, while approximately 10 5 yrs later the gravity of the 
sink particle has altered the gas dynamics thus perturbing 
the shape of the PDF. 

Numerous studies (e.g 



Padoan et al. 



1997 



Passot & 



Vazquez-Semadeni| [1998] |Price et al.||2011| |Molina et al 



2012[ ) have shown the rms Mach number M and gas density 



PDF width a of a supersonically turbulent flow are related 
by 



ln(l + b 2 M 2 ) 



(26) 



where b is a constant, found numerically to lie between ~ 
0.2 — 1.1 (see Federrath et al.|2008[ ). Since ordered rotational 
and radial motions do not directly produce density fiuctu- 
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ations, it is the turbulent Mach number, A4 tU rb = «turb/c s 



with Vturb defined as in Equation ( 24 1, that should be present 
in Equation (|26||. Since the gas is dominated by turbulent 



motions after H2 cooling becomes effective (see Figure|8|, es- 
pecially in low-temperature, self-shielding regions, n rma and 
fturb are very similar here. We can compare the parameters 
of the density PDF of the self-shielding gas in the simulation 
with expectations derived from the independently measured 
turbulent Mach number of the gas flow. Roughly 10 yrs af- 
ter sink particle formation, the width of the density PDF 
is a — 0.93, corresponding to the green curve in Figure |5J 
and the turbulent Mach number in the self-shielding core is 
-Mturb ~ 2.7. Equation (261 would then require b w 0.43, 



in excellent agreement with previous numerical experiments 
of idealized turbulence in which turbulence is driven by a 
combination of solenoidal (divergence-free) and compressive 
(curl-free) modes, i.e., mixed-mode forcing (Federrath et al 
|20T0| ). This seems to suggest that the gravitational infall 
which is driving supersonic turbulent motions in the cold, 
self-shielding core, is exciting a combination of solenoidal 
and compressive modes, and that turbulence has modified 
the density structure of the gas in close accord with theo- 
retical expectations. 

We stress that log-normal density PDFs have been 
found in extremely idealized simulations where supersonic 
turbulence is driven by random forcing. In the simulation 
presented here, however, the turbulence is driven by gravi- 
tational infall and thermal instability. Thus, it is interesting, 
though not surprising, that the parameters of the gas den- 
sity PDF in our cosmological simulation so closely match 
the idealized, theoretical expectations. 



5 GAS FRAGMENTATION AND STAR 
FORMATION RATE 

In this section we discuss trends towards gravitational frag- 
mentation in our simulation. We also attempt to predict 
properties of the expected sub-grid star formation in the sink 
particle, including whether fragmentation on scales smaller 
than the sink particle accretion radius is expected. Our anal- 
ysis is the first step towards learning about the mass spec- 
trum of the stellar objects produced and the overall effi- 
ciency in which this system converts gas into stars. 



5.1 Star Formation Rate and Global 
Fragmentation 

For 3 x 10 5 years after sink particle creation, the extent 
of our simulation, no additional sink particles form. At the 
end of the simulation the lone sink particle has a mass of 
« 1100 M©, which is ~ 6% of the mass of the self-shielding, 
cold gas. To understand this lack of fragmentation and gain 
insight into the expected rate of star formation, we define 
the dimensionless star formation rate per free-fall time as 



SFRff 



M s i n k 

(M/iff) ' 



(27) 



where M and tg are the total mass and characteristic free- 
fall time of the self-shielding region, and Mi n k is the sink 
particle accretion rate (see Figure |4| . SFRff is a measure 
of the actual rate of star formation compared to the rate 



if gas is collapsing on the free-fall timescale. Even though 
M s i n k is the sink accretion rate, not the star formation rate, 
it can be interpreted as the rate at which gas collapses to 
high densities where it can fuel star formation. Given that 
we do not attempt to model any feedback effects from the 
sink particles, such as protostellar outflows and radiation, 
Mi n k is a firm upper limit to the actual star formation rate. 

Taking fiducial values for the self-shielding core, namely 
the average sink particle accretion rate, characteristic self- 
shielding cloud density, and its mass Equation ( 27 1 becomes 



SFRff 



0.1 



M sin 



3 x IO- 3 M Q yr- 



tg 



10 6 yr/ V10 4 M ( 



M 



(28) 



which suggests that star formation in our system is at least 
10 times slower than it would be if all the gas in the self- 
shielding core gas had been collapsing at the free-fall rate. 
We will argue this retardation stems from the supersonic 
turbulence and a centrifugally supported disc around the 
sink particle. 

One way to understand how supersonic turbulence sup- 
presses the rate of gaseous collapse is by considering it an 
additional source of pressure support. The Jeans mass, Mj, 
scales with sound speed and density like Mj oc p _1 ' 2 Cg eff . 
The presence of supersonic turbulent motions enhances the 
effective sound speed, c a off = c 2 + v 2 ms /3. Additionally, 
strong isothermal shocks impart density fluctuations scal- 
ing as p oc M 2 oc t>rms- These two effects lead to a scaling of 
Mj oc v 2 ms in supersonic gas, implying turbulence does in- 
deed impede collapse by increasing the effective Jeans mass 



(e.g., Chandrasekhar| 195 1 Bonazzola et al.|1987 Mac Low 



& Klessen 2004 1. It should be noted, however, that this anal- 



ysis assumes that the highly anisotropic turbulent velocity 
field is isotropic; thus treating supersonic turbulence as an 
additional source of pressure is not necessarily justified. 

It has long been known that the global SFRff in the 
Galactic giant molecular clouds (GMCs) is on the order of 
(e.| 



0.01 



Zuckerman & Evans 19741. This low rate of 



star formation extends to more compact star-forming sys- 
tems as well, such as infrared dark clouds (IRDCs) and 



higher density clumps embedded within (see Krumholz & 



|Tan|2007| and references therein) . Equation ( 28 1 suggests a 
global SFRff of ~ 0.1 for self-shielding, cold gas seen in our 
simulation, ten times that of the Galactic norm. Perhaps 
this is not surprising, though, given we neither include any 
sort of stellar feedback nor magnetic fields, both of which 
would likely reduce the star formation rate. Additionally, the 
rms Mach number in the self-shielding gas here, A4 ~ 3, is 
much less than typical Mach numbers in Galactic molecular 
clouds, M ~ 10 - 20 (e.g., |Bergin fc TafaUa||2007) imply- 
ing a relatively smaller contribution of turbulent pressure 
support in the simulation compared to molecular clouds. 



Indeed, Krumholz & McKee (20051 showed that the den- 



sity threshold for star formation, understood as the density 
where thermal pressure equals turbulent pressure, in a su- 
personically turbulent medium is p CT it ~ po M 2 ■ 
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Figure 10. Cooling timescale t coo i and 3/Q (top panel) and the 
Toomre Q parameter (bottom panel) plotted as a function of 
distance from the sink particle for gas with n > 10 6 cm -3 which 
roughly selects gas in the sink particle's disc. In the top panel, 
thick lines are the cooling time while the thin lines represent 3/Q. 
Both panels suggest that the disc is stable against fragmentation. 
Specifically, in the top panel, 3/C is generally below t coo i (the 
Gammie criterion), while in the bottom panel, Q > 1 throughout 
the disc's extent. For reference, the disc radius is approximately 
~ 0.05 - 0.07 pc. 



5.2 Disc Fragmentation 

The idea that supersonic turbulence suppresses star forma- 
tion and fragmentation does not apply in the ~ 1 pc vicinity 
of the sink particle where the densities exceed ~ 10 6 cm~ 3 . 
In Figure [5] we show that a disc of size ~ 0.1 pc forms 
around the sink particle with rotational velocities approach- 
ing ~ lSkms^" 1 at r ~ 10~ 2 pc. It is interesting why this 
disc does not become unstable and fragment into multiple 
sink particles, given that disc fragmentation has been seen 



in recent simulations of Pop III. 1 star formation (e.g., Stacy 
et al.|2010| |Clark et al.|2011a| [Greif et al.|20fTj ), though at 
much higher densities and smaller spatial scales than probed 
here. 

In general, two conditions must be met for a disc to frag- 
ment. The disc must be gravitationally unstable which re- 
quires that the Toomre Q parameter, Q = CsQ/ttGT,, be less 
than unity ( |Toomre|p6H |Goldreich fc Lynden-Bell||1965| ), 
where c B is the sound speed in the disc, fi = v TOt /R is the 
angular velocity, and £ is the gas surface density. Addition- 
ally, even if a disc is gravitationally unstable, a number of 
effects can prevent fragmentation into gravitationally bound 



collapsing clumps (e.g., Kratter & Murray-Clay 2011). In 
particular, |Gammie| ( |2001[ ) argued that a key effect influ- 
encing the fragmentation behavior of a disc is the balance 
between cooling and the dissipation of turbulence, which is 
a source of heating. Quantitatively, Gammie found that if 
^cooi < 3/0 then discs tend to fragment into clumps. 

In the top panel of Figure[lO]we show the mass- weighted 
average of the cooling time as a function of distance from 
the sink particle for gas with density n > 10 6 cm -3 at three 
different times (thick lines). We also show the mass- weighted 
average of 3/0 (thin lines). As is shown, 3/0 and t coo i trace 
each other in the inner disc (r < 0.03 pc), while 3/0 is a 
factor of ~ 2 lower than t coo i at larger radii and particu- 
larly in the last two times shown. According to the results 



of Gammie ( 2001 1 , this disc is stable to fragmentation given 



its inefficient cooling and relatively rapid rotational velocity. 
For reference, the radius of the disc is sharp and is approx- 
imately ~ 0.05 - 0.07 pc. 

In the bottom panel of Figure[lOj we show the Toomre Q 
parameter as a function of distance from the sink particle. To 
calculate Q(r), we rewrite the gas surface density as E(r) = 
p(r)H, where H is the characteristic vertical thickness of the 
disc. We have verified that for the disc here, H m 0.05 pc is 
approximately constant both in radius and time, and that a 
more detailed calculation of £(r) would come within a factor 
of two of this estimate. Thus, Q(r) provides an additional 
reason for the absence of disc fragmentation since Q > 1 
within the disc's extent and actually increases slightly with 
time. 



6 HD COOLING AND POP III.2 STAR 
FORMATION 

Elevated electron abundances created by the ionizing radia- 
tion emitted by the first stars, or in our case, by collisional 
ionization in TUr > 10 4 K haloes, are thought to result in 



an increase of the HD abundance (see Section 2.3 1. HD is 
an agent that can act as a coolant in metal-free gas below 
200 K. These lower temperatures reduce the Jeans mass and 
may result in stars with lower characteristic masses, the so- 
called Pop III. 2. However, as discussed in Section [3] HD is 
never a significant coolant. Some gas does cool to tempera- 
tures < 100 K as seen in Figure |3j but this is the result of 
adiabatic expansion, not HD cooling. 

In Figure [TT] we show a phase plot of the HD cooling 
time, tcooi, hd as a function of density at the time when the 
sink particle formed. We also show the free-fall time tg. For 
HD cooling to be significant, it is necessary that £ C oo1.hd < 
tff. This criterion is never met. 

Many studies have explored the HD cooling mode in 
Pop III. 2 star formation. Johnson & Bromm (20061 cal- 



culated that an HD abundance of ihd ~ 10 is needed 
for HD emission to cool gas to Tqmb in a Hubble time. 
They suggested that this sets a firm lower limit on the HD 
abundance needed for its cooling to have a thermodynamic 
impact. Using one-zone thermodynamic models, they also 
found that this HD abundance is realized in the downstream 
of T ~ 10 4 K virialization shocks, in supernova remnant 
cooling, and in the collapse of relic Pop III. 1 HII regions. In 



a cosmological setting, Greif et al. ( 2008 \ detected gas with 
temperatures ~ Tomb, which they attribute to efficient HD 
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Figure 11. HD cooling time and the free-fall time as a func- 
tion of density when the sink particle forms. For HD cooling to 
be thermodynamically significant, the HD cooling time must be 
shorter than the free-fall time, which is never realized. This is a re- 
sult of many factors, including the reduction in the free electron 
fraction when H2 cooling becomes significant at n m 200 cm -3 
and the reduction in the HD cooling rate when it thermalizes at 
n ~ 10 5 cm -3 . 



cooling, in an atomic cooling halo of similar mass to the one 



we analyze. Yoshida et al. (20071 found that HD cooling was 



responsible for gas cooling to the CMB temperature after 
the collapse of a Pop III. 1 relic HII region and argued this 
was responsible for reducing the mass of the clump which 
eventually will undergo runaway collapse. There has even 
been a suggestion, by |McGreer fc B ryan ( 2008]) , that in Pop 
III. 1 star formation HD cooling may play a thermodynam- 
ically significant role. In contrast, other studies that have 
examined the collapse of metal-free gas in a similar context 
to the one explored here are in agreement that HD never 
plays a major thermodynamic role | Omukai|2001 O'Shea & 
|Norman||2008| |Shang et al.|2010[ ). 

We propose two related reasons for the lack of sig- 
nificant HD cooling in our simulation. First, as found by 



Wolcott-Green & Haiman (20111 using one-zone models, 



the LW flux necessary to suppress HD cooling is J cr it,HD ~ 
10~ 22 ergs" 1 cm -2 Hz -1 sr _1 . This is far smaller than the 
flux required to completely suppress H2 formation and cool- 



ing, Jcrit.iio ~ 10 



10" 



' erg s cm Hz sr ( Shang 



et al. 20101. This critical HD flux is not the result of di- 



rect HD photodissociation, but is instead due to the partial 
photodissociation of H2 which never permits a significant 
fraction of gas to cool below ~ 300 K where significant HD 
fractionation can occur. Second, at n CO oi ~ 200 cm -3 where 
H2 cooling becomes effective, the free electron abundance is 
not high enough to drive up the H2 abundance. Gas heated 
to high temperatures in virialization shocks can have a free 
electron abundance of x c ~ 10~ 2 , but only at low densi- 
ties, 10~ 2 — 10 _1 cm -3 . At larger densities around n coo i the 
electron abundance drops to ~ 10 -4 (see Figure [7]), below 
the residual abundance after recombination, simply because 
the electron recombination time is inversely proportional to 
density. This would seem to suggest that Pop III. 2 star for- 
mation is very similar to Pop III. 1 , at least in the scenario 
explored here. 



7 DISCUSSION 

7.1 Lyman- Werner Radiation Field 

A key physical parameter in this study is the intensity of 
the constant LW radiation background. Given the large 
mean free path of photons with energies below 13.6 eV, the 
formation of the first stars should have established a LW 
background well before the completion of reionization (e.g., 
Haiman et al.|2000[ |. 

In terms of the cosmic star formation rate, we can esti- 



mate Jlw,2i as (e.g., Johnson||2011 \ 



Jlw,: 



lO -2 M yr -1 Mpc~ 



1 + 2 
10 



,(29) 



where ?7lw is the number of LW photons produced per stellar 
baryon, p\ is the star formation rate per comoving volume, 
and the lifetime of stars producing the LW radiation is im- 
plicitly assumed to be 5 x 10 6 yrs. The star formation rate 
is normalized with a reasonable value around z = 10 (e.g., 
Trenti fc Stiavelli|2009 1. 

Another estimate of Jlw, 21 can be obtained by consider- 
ing the ionizing background needed to reionize the Universe. 
Assuming that reionization was driven by photons with en- 
ergies just above the Lyman limit, it is straightforward to 
estimate the value of t he UV background jus t below the 



limit, J 21 « Jlw, 21 (e.g., Bromm fc Loeb|2003 l, 



400 



N_ 

To 



fese 

oT 



10 



(30) 



where iV 7 is the number of ionizing photons-per-baryon re- 
quired to reionize the Universe and / eS c is the ionizing pho- 
ton escape fraction from high-redshift star forming galaxies. 
Both of the above estimates suggest that before reionization 
was complete, radiation intensities as high as Jlw, 21 ~ 100, 
the value we adopt here, were realized in the Universe. De- 
pending on the scenario for reionization, these may overes- 
timate the mean background intensity, but Dij kstra et al.| 
(2008) showed that, due to dark matter halo clustering, a 
small fraction of haloes can experience LW intensities greatly 
above the global mean value. 

Even though the UV background intensity around the 
time of reionization is very uncertain, the above estimates 
suggest it was likely far above the threshold, Jlw, 21 ~ 10 , 
that would have suppressed the formation of H2 in small- 
mass ~ 10 5 — 10 6 Mq haloes, thus delaying cooling and star 
formation until the assembly of atomic cooling haloes. As 
we shall discuss in Section |7.3| the precise value of Jlw, 21 
may significantly change the results of our simulation, par- 
ticularly the star formation efficiency. 

For more accurate estimates of Jlw, 21 we can look to 
other studies which have undertaken detailed modeling of 
the redshift evolution of the LW background intensity. In 
Figure [L2| we show three such estimates in the literature. We 
also show the LW intensity which would have accompanied 
reionization at z — 9 (top solid line) assuming the fiducial 
values in Equation (301. Additionally, we plot the redshift- 
dependent LW radiation intensity JLW,crit(z) needed to com- 
pletely suppress H2 formation in haloes with T vlr < 10 4 K 
(bottom solid line). We compute JLw.crit(z) by equating the 
JLW,2i-dependent minimum halo mass needed for H2 cooling 
with the halo mass at which the virial temperature equals 
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Figure 12. Theoretical estimates of the redshift evolution of the average radiation intensity in the LW bands (JlW 2l) compiled from 
the literature. The curves show the results of |Ahn et al.| l [2"009| (green triple-dot dashed fmerf), [Trenti fe Stiavelli] ( |2009} (blue dashed-dot 
line), and[Holzbaucr & Furlanctto ( 2 012| ) (red dashed lines). The three curves of|Holzbaucr & Furlanctto (2012J correspond to different 
choices for the minimum virial mass of a halo capable of forming stars, which they take, from top to bottom, as 10 6 , 10 7 , and 10 8 Mq. 
The top black line denotes the LW radiation background which accompanies reionization of the IGM at z = 9, assuming an escape 
fraction of ionizing radiation of / eS c = 0.1 and 10 ionizing photons-pcr-baryon required for reionization. The bottom black line is an 
estimate of the LW background needed to suppress cooling in haloes which rely exclusively on H2. See the text for further details. 



10 4 K (equations 12 and 14 from |Trenti fc Stiavelli|[2009 
respectively) to obtain 



JLW,crit(z) = 1.5 



1 + Z 

31 



(31) 



We refer the reader to the original references for the details 
of the models, but briefly summarize the critical compo- 
nents. Trenti & Stiavelli ( 2009 1 and Holzbauer & Furlan- 



etto 



(20121 model Jlw,2i(z) in a semi-analytic fashion, us- 



ing the formalism of Press & Schechter (1974k modified for 



ellipsoidal collapse by Sheth & Tormen (20021 to estimate 



the mass function of collapsed haloes, which they combine 
with a prescription for star formation in these haloes. |Trenti| 
& Stiavelli (120091) self-consistently utilize their calculated 



LW background to estimate the minimum mass of a star- 
forming halo. The three curves in Figure [12] from Holz bauer] 
& Furlanetto (20121 correspond to different choices for the 
minimum mass of a star-forming halo; M m i n = 10 6 , 10 7 , 
and 10 8 Mq from top to bottom, respectively, each with a 



fixed star formation efficiency of 0.1. Ahn et al. (20091 do 



not rely on a semi-analytic approach, but instead carry out 
transfer of ionizing and LW radiation in a large, cosmologi- 
cal, JV-body simulation. Their resolution, however, is limited 
to 10 8 M©, which may be responsible for their Jlw,21 esti- 
mate being the lowest at high redshifts, as can be seen in 
Figure [12] 

There are numerous uncertainties in these estimates of 
Jlw,2i(z), particularly in the choice of star formation effi- 
ciency (which is likely to depend on Jlw,2i itself, as well as 
the redshift, halo mass, and metallicity), LW escape frac- 
tion of radiation, mass and multiplicity of Pop III sources, 
and different prescriptions for star formation feedback. We 
note that all these models assume that Pop III sources 



have an extremely top-heavy IMF with a characteristic mass 
~ 100 Mq. In spite of these uncertainties, most of these es- 
timates seem to converge around z = 10, but they still un- 
derpredict the LW intensity that would have accompanied 
reionization by z = 9. Additionally, all models suggest that 
at z < 15 — 20, the average LW radiation field is large enough 
to completely suppress cooling and star formation in haloes 

with T vir < 10 4 K. 

Furthermore, as discussed in Ahn et al. ( 2009 1 and 



Holzbauer & Furlanetto (20121 (see also Dijkstra et al. 
2008 1 , the LW radiation field is expected to be very spatially 



uniform around the time of reionization, in the sense that 
isolated regions with Jlw,21 significantly above or below the 
mean were extremely rare. This homogeneity increases with 
decreasing redshift and increasing spatial scale. Compared 
to this spatially homogenous LW radiation feedback, chemi- 
cal enrichment of the IGM by galactic outflows has been sug- 



gested to be extremely inhomogeneous (e.g., Scannapieco et 



al.|2002||Tornatore et al.|2007| [Trenti et al.|2009| . This sug 
gests that the scenario we explore in this work, a chemically 
pristine atomic cooling halo subjected to a strong LW back- 
ground, is a physically plausible scenario for the formation 
of metal-free stellar clusters. 



7.2 Internal Feedback 

This work focused on the impact of an external LW radia- 
tion background and its effect on star formation in metal- 
free atomic cooling haloes. However, once a starburst be- 
gins, LW and photoionizing radiation feedback from inter- 
nal sources will likely impact subsequent star formation. 



Glover & Brand (2001 \ explored the conditions under which 



a metal-free gaseous clump exposed to LW radiation can 
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Figure 13. Gas density as a function of distance from the sink 
particle 10 5 years after sink particle formation colour coded by 
the mass-weighted H2 abundance. The required density for a 
given clump's free-fall time to equal its H2 dissociation time is 
also shown by the solid lines, that were computed assuming all 
clumps have a mass of 3000 Mq RJ Mj and and H2 abundance of 
xn 2 = 10 — 3 . The different solid lines correspond to three different 
choices for the LW photon production rate N^i s : 10 49 s _1 (top) 
is appropriate for a sing le 100 Mq star, lO^s" 1 (middle) for ten 
10 Mq stars, while 10 47 s — 1 (bottom) would be the rate from ten 
5 Mq stars. Assuming a homogeneous clump, its density must lie 
above a solid line for its collapse to occur before its supply of H2 
is exhausted by photodissociation. 



continue to collapse and form stars. They argue that star 
formation in gas clumps will be suppressed if the H2 pho- 
todissociation time, tdis, is shorter than the free-fall time 
of the clump. Glover & Brand (20011 derived that tdis for 



a spherical, homogeneous cloud of density n, mass M, and 
distance D from a source which produces Ndis LW photons 
per second is given by 



t die = 20 a; H 2 n 2/3 M s 1 o / 1 ^^ c / d -: 3 1 / a b 1 s 



( iVdis 

V 10 48 s" 



yrs (32) 



where / a b s is the fraction of incident photons which are ac- 
tually absorbed by the cloud, /dis is the fraction of H2 ex- 
citations by a LW photon which result in a successful dis- 
sociation, and Afsdar and D pc are measured in units of Mq 
and parsecs, respectively. Considering the simulation here, 
a clump at a distance of ~ 1 pc from the source of LW pho- 
tons (which we take to be the sink particle) has a character- 
istic mass on order of the Jeans mass ~ 3000 Mq , density 



10 cm , and an H2 abundance 



10" 



If we reason- 



ably assume that the central sink particle forms roughly ten 
10 Mq stars (see Section [5| the LW photon production rate 
is iVdis ~ 10 4S s _1 . The resulting H2 photodissociation time 
tdis ~ 6 x 10 4 yrs assuming that / d i s = / a b s = 0.1. The 
free-fall time is iff (10 5 cm" 3 ) w 1.5 x 10 5 yrs, therefore this 
clump will have its primary cooling agent, H2, photodisso- 
ciated and thus its further collapse suppressed. 

In reality, the effect of LW feedback on nearby collaps- 
ing clumps will depend strongly on their distances from the 
radiation source and on the level of pre-condensation in the 
gas, i.e., how dense the clump is when the neighboring LW 



source turns on. In Figure[13]we show the radius and density 
dependence of the H2 abundance in cells with ih 2 > 10~ 5 
at a time 10 years after sink particle formation, roughly 
the Kelvin-Helmholtz time of massive stars that would be 
the primary producers of LW photons. The cells are colour 
coded by their H2 abundance. We also show the density 
where iff = tdis for three different choices of the LW pho- 
ton production rate. To derive this density, we assume that 
collapsing clumps have a total mass of 3OOOM Qh 2 abun- 
dance of 10 -3 , and that /dis = /abs = 0.1. Evidently the 
effect of LW radiative feedback on subsequent star forma- 
tion is a strong function of JVdis and the source distance. If 
star formation which occurs in the sink particle results in 
a LW photon production rate of iVdis ~ 10 48 s _1 , it is clear 
from Figure [13] that further collapse and star formation will 
be strongly suppressed in the uncollapsed fraction of the 
self-shielding core, even given the large simplifications we 
have made in this analysis. This may have an adverse effect 
on the overall star formation efficiency and detectability of 
these systems. 



7.3 Direct JWST Observations of a Metal-Free 
Stellar Population 

Since the basic characteristics and formation mechanisms 
of the first cosmic structures are very uncertain, one of the 
main goals of future observations with the JWST is to ob- 
serve light from the first stars and galaxies. By modeling 



metal-free synthetic stellar spectra, Zackrisson et al. (20111 
determined JWST exposure limits and colour criteria for 
high-redshift, metal-free galaxies. Given the uncertainties of 



these sources, IZackrisson et al. (20111 explored a large pa- 



rameter space defined by the shape of the stellar IMF, neb- 
ular emission strength, stellar population age, stellar mass, 
and formation redshift. Their most optimistic scenario in- 
volves a 3 Myr old stellar population with an extremely top- 
heavy IMF (dn/dM oc M -2 ' 35 for 50 < M/Mq < 500) and 
maximal nebular emission, requiring a vanishing escape frac- 
tion of ionizing radiation, / csc = 0. At z = 10, a cluster with 
stellar mass M* « 10 5 Mq and the above characteristics 
would just be detectable (at 10cr) with JWST in ultra-deep 
(100 hr) broadband exposures. 

Given our finding that ss 1.5 X 10 4 Mq of cold gas be- 
comes available for star formation in an atomic cooling halo 
at z ~ 12.1 exposed to a LW intensity of J21 = 100, it 
seems extremely unlikely this starburst could be detected 
by JWST, even if 100% of the mass turns into stars with an 
extremely top-heavy IMF. It is conceivable, however, that 
a weaker UV background could still have suppressed star 
formation until the assembly of an atomic cooling halo, but 
could have allowed for a larger global star formation effi- 
ciency. An inspection of Figure [12] suggests that a back- 
ground LW intensity Jlw,21 > 1, at least for z < 30, would 
have fully suppressed star formation in halos not capable of 
atomic cooling. 



We note that the clumps visible in Figure |13| at ~ 0.5 pc and 
~ 1.5 pc are in fact not Jeans unstable and have masses well 
below 3000 Mq . This downward mass revision will decrease their 
H2 photodissociation time which depends weakly on mass, t^m oc 
MVS. 
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To get a sense of how the mass of cold gas depends on 



J21, we proceed as follows. First, recall Equation ( 18 1 which 
expresses the density at which H2 cooling becomes effec- 
tive, n coo i, as a function of J21 and x e . This cooling den- 
sity partially determines the mass of cold T < 10 3 K gas, 
A/coid, available for the first burst of star formation. Once 
a halo reaches the atomic cooling limit, the accretion rate 
of gas onto the self-shielding core, Af co id, will not a strong 
function of J21. With this, we can estimate the cold gas 
mass as M co id ~ Mcold te(^cooi) oc J^j 1/ ' 3 xl . The abun- 
dance of free electrons would be higher at lower densities, 
but by how much is difficult to determine given the highly 
non-equilibrium electron abundance. For example, if x c is 10 
times higher than in our simulation at n coo i and if J21 = 10, 
the mass of cold gas would be M coM ~ 10 2/3 x 1.5 x 10 4 Af© « 
7x 1O 4 M . If we furthermore assume a relatively high star 
formation efficiency, /„ = 0.1, which represents the fraction 
of cold gas that will eventually turn into stars, we find a to- 
tal stellar mass of A/* ~ /* M co id ~ 7000 Af . Even if these 
stars form with an extremely top-heavy IMF, a possibility 
recent simulations suggest is unlikely (e.g., |Stacy et al.|2010| 
Clark et al.|2011a Greif et al.|2011 1, the analysis of Zackris- 



son et al. (20111 implies this cluster will not be detectable 



with the JWST. 

One effect, however, which may render these targets 
within the JWST detection limit is gravitational lensing. 
Low-redshift galaxy clusters are capable of magnifying high- 
redshift sources by factors of fj, ~ 10 — 100. Indeed, the re- 
cently disc overed z « 9.6 galaxy MACS 1149-JD1 ( |Zheng 
et al. 2012) is gravitationally lensed by a foreground galaxy 
cluster with a magnification factor of /t ~ 14. though lensing 



was not strictly required for its detection. To this end, |Za-| 
|ckrisson et al.| ( |2012[ ) estimated the conditions under which 
the JWST could detect gravitationally lensed, metal-free 
stellar populations behind the z w 0.546 galaxy cluster 
MACS J0717.5+3745, an ideal gravitational lens for high- 
redshift stellar sources. Mock halo catalogs of metal-free 
galaxies, generated with the methodology of |Trenti et al.| 
(2009), were simulated to lie beyond this particular galaxy 
cluster and the expected number of metal-free stellar pop- 
ulations as a function of limiting magnitude was gener- 
ated. The primary source of uncertainty in these estimates 
is e, the efficiency of Pop III. 2 star formation, defined as 
M* lP o P in.2 = eAfhaio(f2b/fi m ). This efficiency could be bro- 
ken down as e = /»/ co id, where / co id is the fraction of a 
halo's baryonic mass which is cold and available for star for- 
mation and /„ is, as before, the fraction of cold gas that will 
eventually become incorporated into stars. 



Zackrisson et al. (20121 find that with efficiencies as 
d , clusters of metal-free stars with moder- 
10 Mq, will be visible in 



low as e ~ 10 

ately top-heavy IMFs, M c h ar 
deep JWST lensing surveys around MACS J0717.5+3745. 
Our optimistic J21 = 10 scenario, described above, results in 
an efficiency e « (7000 M /3 x 10 7 M e )(Q m /Q b ) « 10 -2 9 . 



Thus, according to the estimates of Zackrisson et al. (20121 



and our relatively crude estimate for the Pop III. 2 star for- 
mation efficiency, it appears that JWST will marginally de- 
tect metal-free stellar populations in deep, foreground clus- 
ter lensing enhanced surveys. This conclusion, though, is 
very approximate and more detailed studies are needed to 
fully assess the observational prospects for detecting these 
primitive, low luminosity stellar clusters. 



8 SUMMARY AND CONCLUSIONS 

We have performed a high-resolution cosmological simula- 
tion focused on the collapse of metal-free gas in a region 
of the Universe subjected to a strong molecule-dissociating 
Lyman-Werner background. We resolved densities up to 
10 s cm" 3 and length scales down to ~ 1000 AU. This simu- 
lation utilized a fully non-equilibrium primordial chemistry 
network coupled with a non-local column density calcula- 
tion to accurately compute the H2 and HD photodissocia- 
tion rates. Additionally, our use of sink particles allowed us 
to evolve the simulation for many free-fall times past the 
first gravitational collapse. 

Here we summarize the main findings of this work: 

• With a LW intensity J21 = 100 in a 1 Mpc 3 (comoving) 
box, effective H2 cooling first occurs in a 3 x 10 7 Mq halo at 
z « 12.1. This cooling results in a thermal instability and a 
cold, dense core develops in which H2 shields itself from LW 
radiation. Upon sink particle creation at n = 10 s cm -3 , the 
central core containing self-shielding, cold gas has a mass of 
~ 10 4 Mq, characteristic size ~ 5pc, temperature ~ 400 K, 
mean rms velocity u rms « 6kms _1 , and Mach number w 2. 
Mach numbers increase to « 2 — 4 within 3.5 x 10 s years. 

• Within the inner ~ 10 pc of the halo, the gas flow be- 
comes supersonic only when H2 cooling becomes effective. 
The gas density PDF in the cold, self-shielding gas is approx- 
imately log-normal and has a width consistent with the tur- 
bulent Mach number of the flow given expectations derived 
from idealized simulations in which turbulence is driven by 
a mixture of solenoidal and compressive forcing. With time, 
the density PDF acquires a high density power-law tail as 
self-gravitating gas decouples from the turbulent flow. 

• The rate of gaseous collapse on scales comparable to 
the size of the self-shielding region is suppressed compared 
to what would be expected if all the gas in the self-shielding 
core was collapsing at its free-fall rate. By the time the sink 
particle grows to ~ 6% of the total mass of the self-shielding 
cloud, no other sinks have formed, suggesting any additional 
sites of fragmentation remain unresolved in the simulation. 
We find an upper-limit to the star formation rate per free- 
fall time of SFRff « 0.1 on scales of ~ 10 pc. We argue this 
is due to the additional effective pressure provided by infall- 
driven supersonic turbulence on large scales and centrifugal 
support on smaller ones. 

• HD cooling was found to be thermodynamically in- 
significant for the entirety of the gas evolution. We attribute 
this to partial photodissociation of H2 and the low free elec- 
tron fraction when H2 cooling becomes significant. Given the 
relatively low LW intensity needed for HD formation to be 
suppressed, Jlw,21 ~ 0.1, in comparison with the relatively 
high intensity required to delay gaseous collapse until the 
formation of atomic cooling haloes, Jlw,21 ~ 10, it seems 
unlikely that HD cooling is ever significant in Pop III. 2 star 
formation, at least in the scenario explored here. 

• Given the similarities between the thermodynamic be- 
havior of Pop III forming halos with and without LW feed- 
back, and the complete absence of significant HD cooling in 
both cases, we suggest that Pop III. 2 star formation delayed 
by LW radiation is very similar to Pop III. 1 star forma- 
tion. Additionally, the similarities of our results with other 
simulations that probed higher densities and smaller spatial 
scales (e.g., Stacy et al.|2010 Clark et al.|2011a Greif et al 
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2011) seem to suggest that Pop III. 2 stellar masses in the 



LW-delayed mode are likely comparable to Pop III. 1 masses, 
~ 1 — 40 Mq, possibly even lower. 
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